; docformat = 'rst'
;
; NAME:
;   MR_AreaGrid
;
;+
;
; :Purpose:
;   Calculates grid-cell areas (m^2), given longitude and latitudes of cell center-points.
;   Only for regular grids at present.

; :Categories:
;	  Model utilities
;
; :Returns:
;   2D array of grid cell areas in m^2. Dimensions of n_elements(lon) x n_elements(lat)
;
; :Params:
;	  Lon : longitude of center for each grid cell in degrees
;	  Lat : latitude of center for each grid cell in degrees
;
; :Usage:
;	  Regular grids only.
;
; :Examples:
;	  IDL> area=mr_AreaGrid(lon, lat)
;
; :History:
; 	Written by: Matt Rigby, MIT, May 29th 2011
;
;-
function mr_AreaGrid, lon, lat

  compile_opt idl2, hidden
  on_error, 2

  re=6367500.0	;radius of Earth in m
  
  dLon=abs(mean(Lon[1:n_elements(Lon)-1] - Lon[0:n_elements(Lon)-2]))*!pi/180.d
  dLat=abs(mean(Lat[1:n_elements(Lat)-1] - Lat[0:n_elements(Lat)-2]))*!pi/180.d
  Theta=double(!pi*(90.-Lat)/180.)
  
  AreaGrid=dblarr(n_elements(Lon), n_elements(Lat))
  
  For LatI=0L, n_elements(Lat)-1 do begin
  For LonI=0L, n_elements(Lon)-1 do begin
    if Theta[LatI] eq 0. or theta[LatI] eq !pi then begin
      Lat1=Theta[LatI] - dLat/2.
      Lat2=Theta[LatI] + dLat/2.
      AreaGrid[LonI, LatI]=(re^2.)*abs(cos(dlat/2.)-cos(0.))*dLon
    endif else begin
      Lat1=Theta[LatI] - dLat/2.
      Lat2=Theta[LatI] + dLat/2.
      AreaGrid[LonI, LatI]=((re^2.)*(cos(lat1)-cos(lat2))*dLon)
    endelse
  endfor
  endfor

  return, areaGrid

end